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Abstract 

In this paper we develop a geometric discretization of the EPDiff equa¬ 
tions in one-dimensional case. We extend the method presented in | [2Q| to 
apply to all (not only divergence-free) vector fields and use a pseudospectral 
representation of a vector field. This method can be extended to a multidi¬ 
mensional case in a straightforward way. 


1 Introduction 


The main objective of this paper is to develop a general method of geometric dis¬ 
cretization for infinite-dimensional systems and apply this method to the EPDiff 
equation. Geometric integration has been a very large and active area of research 


(see p8| for an overview). Unlike conventional numerical schemes, geometric in¬ 
tegrators are derived from variational principles and preserve the structure of the 
original systems. The structure-preserving nature of these methods allows to cap¬ 
ture dynamics without usual numerical artifacts such as energy or momenta loss. 

To construct a variational integrator for an infinite-dimensional system, such 
as the EPDiff or Euler equations, one first has to develop a method of discretizing 
the configuration space of this system, i.e. the group of diffeomorphisms. More¬ 
over, we have to replace this group with a finite-dimensional Lie group in order to 
preserve the symmetries of the original system. As the second step we can derive 
a finite-dimensional system on this group from Lagrange-D’Alembert principle. 
Lastly, we apply standard techniques of variational integration to discretize time 
and get an update rule. 

The method described below extends one developed in pO] for incompressible 
Euler fluids. Here this method is presented in a general case applicable to all, not 
only divergence-free, vector fields. Also, a different (pseudospectral) representa¬ 
tion of the velocity field is used. We will apply this method to the one-dimensional 
EPDiff equation and present numerical results in Section]^ 
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1.1 The EPDiff equations 


The EPDiff equations comprise a family of geodesic equations on the group of 
diffeomorphisms Diff(M) of a manifold M, dimM = n, where the metric is 
defined by a norm on the space of vector fields Vect(M) of the following form: 


|2 

\l 


V r. = 


/ {Lv^v)dx. 

Jm 


( 1 ) 


Here (•, •) is the inner product on and L is a positive definite self-adjoint differ¬ 
ential operator. This equation plays a central role in computational anatomy, where 
the distance between an image and a template is measured as a length of a geodesic 
connecting them. See pT| for details. 

Later in this paper we will use iht flat operator instead of L: 

\> \ V ^ ^ {v^^u) = (L'u, u)^ for any v^u ^ Vect(M), (2) 


where Q{M) is the space of one-forms on M and (•, •) is the pairing of a one-form 
and a vector field. 

The EPDiff equations can be derived from the following variational principle: 

s[ [ {v\v)dxdt = 0, = ^1^^^ = ^1^^^ = 0. (3) 

Jo J M 

The constraints on Sv are called Lin constraints in fn\ and are due to the fact that 
the variations are taken along a path on the Lie group Diff (M) while v belongs to 
its Lie algebra. Substituting the expression for 6v into the integral and using the 
fact that the commutator of vector fields is the Lie derivative [v^u] = LyU, we get 



0 JM 


i) + Ly^)^ dxdt = 0, 


which after integration by parts becomes 



0 Jm 


{—v^ — LyV^ — div V, ^)dxdt — 0. 


Thus, we obtain the EPDiff equation: 


& + LyV^ + div 'y = 0. 


( 4 ) 


Later on in this paper we will consider a special case of the EPDiff equation 
when dim M = 1 and Lv — v — this case the EPDiff equation becomes 

the Camassa-Holm (CH) equation: 

m+{mv)x + mvx = 0, m = v - ( 5 ) 
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which is a well known model for waves in shallow water (see Q). This equation 
is completely integrable and has soliton solutions called peakons which have a dis¬ 
continuity in the first derivative. Due to this, solving the CH equation numerically 
can be challenging. 

1.2 Overview of the method 

To construct a discrete version of the EPDiff equation, we will use the method in¬ 
troduced in p0| to discretize the Euler equation of ideal incompressible fluid. In 
this paper, however, we extend this method to apply to the whole space of diffeo- 
morphisms in a pseudospectral representation of the velocity. 

According to this method we replace the group of diffeomorphisms with a 
group of matrices, on which we will construct a Lagrangian system with nonholo- 
nomic constraints. The derivation of the finite-dimensional version of the EPDiff 
equation on a matrix group will closely follow the derivation of the EPDiff equation 
presented above. 

2 General method 

2.1 Discrete diffeomorphisms 

Following po) we will replace a diffeomorphism g € Diff (M) by a linear operator 

Ug: 

Ug : L2{M) ^ L2{M), Ug:((>^ct>og-\ 

where L 2 {M) denotes the space of square-integrable functions on M. We will 
consider a finite-dimensional linear operator q as an approximation to the diffeo¬ 
morphism g and write q ^ g if q approximates Ug. 

To discretize the linear operator Ug we first need to discretize the space where 
it acts, i.e. the space of L 2 functions on M. To do this we fix a family finite¬ 
dimensional spaces C L 2 (M), dim^Ar = N and two families of operators 

Dn • L2{M) Rat • ^N- 

We will call the family Dn a discretization of L 2 {M) if for any function (j) G 
L 2 {M) the sequence (j)N — Daa 0 converges to 0 as W ^ oc. We will call the 
W-dimensional vector = R^^^at a discrete function and the operator Hn a 
reconstruction operator. 

Now we can define a discrete diffemorphism as a linear operator acting on 
discrete functions: 
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Diff(M) 


C^{M) 


C^{M) 



Figure 1: Discretization and reconstruction operators. Here 
is a space of discrete functions and Rat is a bijection. V{M) 
is the group of discrete diffeomorphisms, which is a finite¬ 
dimensional group of linear operators. 


Definition 1. Let Dn be a discretization of L 2 and Ra^ a family of reconstruction 
operators. We will say that a family of linear operators qn : ^ is an 

approximation to a diffeomorphism g G Diff(M) and write Qn 9 if for any 
function f G L 2 {M) we have: 

Raa^atR^^Dat^ ^ Ugfi when N ^ 00 . (6) 

Thus, to discretize the group of diffeomorphisms we first need to choose a 
discretization of L 2 functions and then fix a group of linear operators acting on 
the discrete functions. Different methods can be used for both of these steps, we 
will describe one such method in more detail below. After the set of discrete dif¬ 
feomorphisms has been chosen we will denote it V{M). The relationship between 
V{M) and Diff (M) is illustrated by the diagram in Figure[^ Note that the diagram 
doesn’t commute. 

2.2 Discrete vector fields 

To define a discrete vector field let’s consider a smooth path qt G V{M) of discrete 
diffeomorphisms. A discrete function is transported by the flow qp. 

It satisfies the equation 

= Qt4>i = = Utp, (7) 

where Ut = Note, that this equation is analogous to the advection equation 

ft — lit 0^5 
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where is the Lie derivative along the vector field Uf. Thus, the linear operator 
Ut — can be considered a discretization of the Lie derivative, which brings 

us to the following definition: 

Definition 2. Let be a discretization of L 2 and a family of reconstruction 
operators. We will say that a family of linear operators Un : ^ is 

an approximation to a vector field u G Vect(M) and write Un u if for any 
function 0 G (7^ (M) we have: 

^ —Luf^ when N ^ 00 ^ (8) 

where convergence is assumed to be in L 2 norm. 

Now, if we assume that the discrete diffeomorphisms V from a Lie group, we 
can see that the space of discrete vector fields, which we will denote by 2), is the 
Lie algebra of V. Moreover, the conunutator [U^V] = UV — VU of two discrete 
vector fields is an approximation to the commutator of the continuous vector fields 
u and V, assuming U u and V ^ v. If the space of discrete functions ^ 
has dimension N, the space of discrete vector fields may have dimension as large 
as N‘^. To make the discretization computationally tractable we will restrict the 
discrete vector fields to belong to a space S of dimension 0{N) instead. However, 
the space S is likely not closed under commutators, [5, S] ^ 5, and therefore we 
cannot restrict discrete diffeomorphisms to a subgroup of V. A method to construct 
a constrained set S will be outlined below. 

For every vector field v G Vect(M) we will be able to construct its discrete 
version 1/ G 5, thus we will define an operator S : Vect(M) ^ S. We will require 
this operator to be right-invertible, so any matrix V ^ S can be reconstructed 
into a vector field. Later in this paper we will use a pseudospectral representation 
in which a vector field on a circle is represented by its values at N points. The 
operator S will be defined in ( [3T] ). 

Note that the matrices in the commutator space [5,5], however, cannot be 
identified with continuous vector fields. See figure]^ 

2.3 Discrete forms and flat operator 

Let’s assume the space Vect(M) is equipped with an inner product (•,•). A discrete 
version of this inner product can be defined as follows: 

Definition 3. A family of Hermitian forms on T)n is said to be an approx¬ 

imation to the inner product (•, •) if for any pair of vector fields u, v ^ Vect(M) 
and its discretization Un u, Vn ^ v, such that Un ^ S, Vn G 5 U [5,5] we 
have 

{Un’>Vn)% ^ when 00 . (9) 
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C'°°(M) > C^{M) 



Figure 2: Discretization of vector fields. Here 5 is a set of lin¬ 
ear operators on representing vector fields and D is the Lie 
algebra of the group of discrete diffeomorphisms V. 


Later on we will omit the superscript d in the formula above and simply write 
(U^V) for the discrete inner product. 

An inner product (•, •) on Vect(M) defines aflat operator 

\) :u^ vt ^ Q}{M)^ (u^v) = for any 'L’ G Vect(M), 

where is the space of one-forms on M. 

Following p0| we define a discrete one-form as an object dual to the discrete 
vector fields, i.e. as a matrix F and a pairing 

{F,U) =Tr(Ff/*). 

This definition of the pairing allows us to define a discrete flat operator \) \ U 
as 

\):U (U, V) = Tr(C/V*), for any ^ G Vect(M). (10) 

2.4 Lagrangian mechanics on the group of discrete diffeomorphisms 

Our goal is to construct a Lagrangian system on the group V{M) of discrete dif¬ 
feomorphisms approximating a certain continuous dynamics on Diff(M). To do 
this, we will construct a Lagrangian of the form (see section |3.2| for an explicit 
construction of the flat operator) 

L{U) = \{U\U) (11) 

and derive the dynanfics from the Lagrange-D’Alembert principle: 

5 [ L{U)dt = 0, 5qq-^ E S,U e S, 5q(0) = 5g(l) = 0. (12) 

Jo 
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The equations describing the dynamics can be easily derived as follows: first, since 
U — qq~^ we can show that 6U has to satisfy the Lin constraint: 

6U = S + [f/, S], where B — 5qq~^. (13) 

Second, substituting the Lin constraint into the expression for 5L{U) we get 

5L{U)^]^{U\B+[U,B]). (14) 

Thus the Lagrange-D’Alembert principle may be written as 

[ Tr{U^*{B+[U,B]))dt = 0, for any B G S', B\^^^ = B^^^ = 0, 

Jo 

which after integration by parts and rearrangement by permuting under the trace 
yields 

+ [U*,U^],B) = 0, for my BeS. (15) 

2.5 Discrete time 

To discretize time we consider the dynamics is given as a discrete path qo^... ^qx 
on V{M), where motion is sampled at regular time intervals tk = k • dt, where 
dt is a time step. For a given pair of configurations qk^ qk+i we use one of the 
following ways to define matrix U for discrete time: 



Qk+l - Qk 

= dt Uk qk, 

(explicit Euler), 


Qk+l - Qk 

— dt Uk qk+i, 

(implicit Euler), 


Qk+i - qk 

= dt Uk 

(midpoint rule). 

(qk+i - q 

,C+1 + C^ 

'fc) 2 

= dtUk, 

(average explicit-implicit). 


These four approaches to discretization result in the following four representations 
of the discretized variational relations: 

1. Explicit Euler. In this case, Uk — {qk^i ~ Qk)/dt q^^. The variation 6kUk 
and 6k^iUk with respect to qk and qk^i respectively become: 

SkUk = 

^k+iUk = —Sqk+iqp. 
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If we denote, similarly to the continuous case, = SqkQj^ we get: 

hUk = + UkBk 

and 

^k+lUk = + ^k+lUk- 


2. Implicit Euler. In this case Uk = ^k+i' H yields: 

^kUk = —-j^SqkQj^^i 

and 

4+lt/fc = ISqk+iQkh - ^^ C+l4fc+lC+l- 

Similarly to the previous case we now obtain: 

hUk = - BkUk, 


and 


h+iUk 


Bk+i 

dt 


UkBk^i^ 


3. Midpoint. The Eulerian velocity between qk and qk+i is now expressed as 

= 2 gfc+i-gfc ■pjjus^ 

“ ((Zfc+i + %) + %) ^ 

~ dt “1“ Qk} 


4. Average Explicit-Implicit. Here the velocity between qk and is ex¬ 
pressed as an average of the velocities computed with explicit and implicit 
rules: ^ ^ 

Uk^ ^-j^{qk+i-qk){qk^+ qkli)- ( 16 ) 

In this case the variations Sk^k+iUk are also averages of the corresponding 
variations: 

SkUk = -^ + l[Uk,Bk], (17) 

^k+iUk+ ■^[Bk+i,Uk]- (18) 
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Now that we have these four different ways to compute variations of Uk, we 
can proceed to derive the corresponding discrete Lagrange-D’Alembert equations, 
e define the discrete-space/discrete-time Lagrangian Ld{qk^ Qk+i) as 

Ld{qk, Qk+i) = L{Uk)> 

The discrete action Ad along a discrete path is then simply the sum of all pairwise 
discrete Lagrangians: 


K-l 

^d{qo, • • •, Qk) = ^ Ld{qk, qk+i)> 

k=0 

We can now use the Lagrange-d’Alembert principle that states that 6Ad = 0 for all 
variations of the qk (for = — 1, with qo and qx being fixed) in Sq while 

Ak is restricted to S. 

Setting the variations of Ad with respect to Sqk to zero for fc G [1,K — 1] 
yields: 

5,(ul_„U,A + ^k{ulUk)^0. (19) 

Now, let’s solve it for Uk in the explicit case. Substituting the expressions for 
SkUk and 5kUk-i yields: 

TT[-UliB*k + dtB*kU*k) + uLi{B*k + dtUUB*k)] = 0 . 

Denoting Ul = {Ul — Ul_^)dt~^ we can rewrite the last equation as 

Tv[{Ut - UlU*k + UUUl-dBl] = 0. (20) 

Let’s fix a basis Bk of the space S, i.e. any matrix [/ G 5 can be written as 

U^Y.^kBk- ( 21 ) 

k 

Now let’s rewrite the equation ( [^ in the coordinates X. First, we have 

rr* TT^ TT^TT^ _ \ ^ '^k—l D* -\rk—l td\> vk td\> vk D* 

Now, if we denote hy A • B the Frobenius product of A and B, we can write 

{uUuLi - ulu^) Xt"X^-\B*B^^) • B, - X^X^iB^B]) ■ B,. 

Let’s denote 

(B;s5) • X = Tr(S*s5s;) = {BiX) • ^5 = (B), BiBp) = 
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and 


■ Bp = TriB^Bp;) = {BpB,) ■ B\ = {B^B^B,) = A,p 

~ ^ip- 

Then the update rule for the explicit case can be written as 

i i 

E CijpXt^X’;-^ - E Di.pXfx’; = 0 (22) 

Similarly, in the implicit case we get 
J2E^pX^ -J2EipX^-^+ 

i i 

E CijpX^X^ - E DijpXt^X^-^ = 0. (23) 

hj hj 

In the average explicit-implicit case the update rule is the average of the two formu¬ 
las above. The midpoint case yields third order terms in U and it’s not considered 
here. 

3 Pseudospectral discretization 

3.1 Discrete functions and vector fields 

To illustrate the method, we consider the following case of pseudospectral dis¬ 
cretization. Let’s define the space S of discrete vector fields on using a pseu¬ 
dospectral representation. Note, that a matrix U G 5 is an approximation to an 
operator of Lie derivative Lu'. 

— 4>cU‘ 

Now we will consider a continuous test function being represented by its trun¬ 
cated Fourier series, i.e. by a vector • • •, where 

4 = ^ f e-^^^cp%x)dx, k = -N,...,N. 

We will denote by D the operator of differentiation in the truncated Fourier space, 

i.e. 

{D4>)k = Dkk^k = ik(t)k- (24) 
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If we know values of a vector field i/(x) at points = —7r+A;-27r/(2A^+l) 
we can define a discrete version of the multiplication operator (j)^ ^ (j)^ • u 

M = FTuF-^, (25) 

where F is the discrete Fourier transform and {Tu)ij = SijUi. 

Now, the space 5 of discrete vector fields is spanned by matrices BkS of the 
form 

Bk = MkD, (26) 

where 

Mk = FIkF~^, {h)ij = ^ij^ik- (27) 

To summarize, our discretization consists of the following: 

1. Space of functions 

■ [(t){x)\^{x)= (28) 

I k=-N ) 

2. Discretization operator: 

Djv : 4>{x) ^ (</>-Jv, • • •, 4>n), ^ y (29) 

3. Reconstruction operator: 

N 

Rat : (0-Ar,..., ^at) ^ 0(x) = ^ (30) 

k=-N 

4. Discretization of a vector field: 

S:vix)^FnF-^D = J2XkBk, (31) 

where 

/u{x-j^) ... 0 ... \ 

Ty = I 0 ... u{xo) ... 0 j , (32) 

\ 0 ... 0 ... u{x]^)j 

Xk^u{xk). (33) 
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3.2 Discrete flat operator 

Let’s now define a flat operator, which is the key ingredient of the method. To 
define a pairing between discrete vector fields U and V let’s note that since U<^ ^ 
Lucj)^ we have for = (0,..., 1,..., 0): 

Uck « = Dkkue^'^^, (34) 

where is defined in the sense of L 2 norm. If a function is represented by a 
vector (j) then 00 ~ / 0- Thus, {Uek)o is m approximation to the — fc-th Fourier 
coefficient of u multiplied by 

{Uek)o = Uok « Dkk j (35) 

Therefore, we can define a flat operator through the following pairing: 

+ = E ^okVok (« - . (36) 

^ -l^kk J-^kk ^ ^ 

It’s worth noting that the pseudospectral discretization allows us to construct a 
flat operator in a much more straightforward way than, for example, discretization 
described in [ [20| . 

3.3 Update rule 

Now, let’s compute the update rule for the explicit case. 

Theorem 1. The update rule in the explicit and implicit cases are given by the 
formulas ( [22| ) and ( [2^ , where 


= ^iF-^DFX)p{F-^HFX)p « 

hj 

1 

(37) 

J^XiXjDijp = ^{F-WF)iX*F-^HFX) « 


(38) 

Y.X,E., = ^(F-'HFX), « im. 


(39) 


i 
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Proof. We have 


Y,XiX^{B),BiBp) = Y,XiX^{B],FTiF-^DFTpF-^D), (40) 

hJ hj 

where = u{xk) (we will write X ^ uin this case). Also, 


(Sj, U) = 5](FI,F-il))o.C7o«(a - D-^) = FojF-^DU<^ - 

S S 

Thus, we can write 


D-")Uos. 

( 41 ) 


5 ] XiXjCijp = 5 ] XiXj {B^j,BiBp) = 

hj hj 

Y,XiXj{FTiF-^DFTpF-'^D)osFojF-\aDss - D-^) = 

^ XiXjFoiF~^DkkFkpFpJ-DssFojFrJ^{aDss - D~J-). 

i,j,s,k 

Since F is unitary and Foi = we have 

F^.fmFu^ = -^(FXU 

i ^ 

and 

Y:x,F,jF-,' = ^(FX).. 

3 ^ 

Now we have 

y] XiXj{B],BiBp) = ^{F-^DFX)p{F-^HFX)p, (42) 

hj 

where 

if = Id - aD‘^. (43) 

Since X u we have 

F-^DFX « 

and _ 

F~^HFX ^ u — auxx- 
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Therefore, now we have 


XiXjCijp ~ ^UxTfl^ 
hj 


(44) 


where m — u — auxx- 

Similarly, 

= Y,X,X,{BlBpB,) = 
hj hj 

J2XiXj{FlpF-^DFljF-^D)osFoiFr/{aD,, - D"/) = 

XjXj(FopF^p DkkFkjFjs ^ 55 )^ 02 ^^5 {o^Dss — ). 

We have 

Y^XiF-'FoiFr'H,, = 


hj j:S 

^{F-'^DF){X * F-^HFX) = ^{F-'^DF){X * F-^HFX), 

where {X ★ Y)i = XiYi. Again, since X « it we have 

Xi.F-'^HFX F^um 


Therefore, 

^XiXjDijp « —dx{um) = — (na;m + Mma;). 


Finally, we compute J2i ^iEip = Xi{B\, Bp): 


(45) 


Y^MBlBp) = - £>-") = 

z i^s 

Y,XiFoiF-^DssFopF-,^DUc^ - D-^) = ^{F-^HFX)p « ^m. 

i^s 


□ 
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Figure 3: Formation of peakons from a gaussian initial condition 


4 Results 

We have implemented our method for the explicit, implicit and the average cases. 
In all numerical tests we see the energy decreasing in the explicit case and increas¬ 
ing in the implicit case. In the average explicit-implicit case however the energy 
is stable. That is, the energy is oscillating around its correct value (see figure]^. 
This behavior is different from energy behavior of a variational integrator. This 
difference is a result of imposing nonholonomic constraints. The same behavior 
has also been observed in other systems of the same form, i.e. for the equation ( [TS] ) 
with a different flat operator. 

We studied different cases of peakon dynamics, such as formation of peakons 
from a gaussian initial condition, interaction of peakons of the same sign and 
peakon collisions. Formation of peakons from a gaussian initial condition is shown 
in Figure For this case we chose a = 1, W = 1000 and dt — 0.01. Peakon 
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Figure 4: Energy behavior for peakon formation 


collision remains a challenge. The simulation leads to creation of multiple peaks, 
but remains stable (see figure]^. The energy drops when the two peakons collide 
initially, but then recovers and remains stable (see figure [^. 



Figure 5: Accuracy and convergence for a single peakon 


5 Conclusions and summary 

To summarize, we have developed a method of discretization for systems on the 
group of diffeomorphisms. This method is presented here for the case of the 
Camassa-Holm equation, but can easily be applied to other systems. The method 
itself is flexible and can use different representations of vector fields (operator S 
in Fig. [^. The final update rule is derived from a variational principle with non- 
holonomic constraints and the resulting energy behavior is different from that of 
a variational integrator. Namely, the energy behavior depends on how the discrete 
velocity is computed from a pair of configurations (see Section [Z5] ). In the average 
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m/m 



Figure 6: Energy behavior for three different methods 


explicit-implicit case the energy remains stable over long time. 

6 Future work 

While the time-continuous system ( fT5] ) is energy-preserving, the energy behavior 
of the time discrete system depends on the choice of discretization of U. One may 
use an adaptive time step method described in Q to construct an energy-preserving 
integrator. However, the effect nonholonomic constraints have on a variational 
integrator remains an open question. 
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Figure 7: Peakon collision sequence shows unstable behavior 



Figure 8: Energy behavior for peakon collision. Energy jumps 
during the initial collision but returns to the neighborhood of its 
correct value after that. 
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